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Abstract 

In order to precisely quantify the fundamental interactions between heterogeneous 
lipid membranes with coexisting liquid-ordered (Lo) and liquid-disordered (Ld) do¬ 
mains, we performed detailed osmotic stress SAXS experiments by exploiting the 
domain alignment in raft-mimicking lipid multibilayers. Performing a Monte Carlo 
(MC) based analysis allowed us to determine with high reliability the magnitude 
and functional dependence of interdomain forces concurrently with the bending 
elasticity moduli. In contrast to previous methodologies, this approach enabled 
us to consider the entropic undulation repulsions on a fundamental level, without 
having to take recourse to crudely justified mean-field like additivity assumptions. 
Our detailed Hamaker coefficient calculations indicated only small differences in 
the van der Waals attractions of coexisting Lo and Ld phases. In contrast, the re¬ 
pulsive hydration and undulation interactions differed significantly, with the latter 
dominating the overall repulsions in the Ld phase. Therefore, alignment of like 
domains in multibilayers appears to originate from both, hydration and undulation 
repulsions. 

Key words: osmotic stress experiments; interbilayer forces; liquid ordered phase; 
liquid disordered phase; membrane rafts; Monte Carlo simulations 
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Introduction 

Diverse physiological processes in living systems depend on fundamental physical 
interactions between lipid membranes acting on the nanoscopic length scale. Of 
particular interest in this context are, besides intra-membrane interactions 
forces acting between membrane domains/rafts across the aqueous phase, which 
are also involved in their correlated mutual alignment. Such positional correla¬ 
tions are well established for liquid-ordered (Lo)/liquid-disordered (Ld) domains 
in model lipid multibilayers (□hud- Several groups have established composi¬ 
tional phase diagrams for mixtures of high-melting lipid, low-melting lipid and 
cholesterol, which exhibit Lo/Ld phase coexistence over a broad range of compo¬ 
sitions and temperatures dUd. These systems mimic mammalian outer plasma 
membranes and enable studies of domain properties under well-defined conditions. 
Most recently, we reported structural details of Lo/Ld phases in two ternary lipid 
mixtures using a global small-angle X-ray scattering (SAXS) analysis for coexisting 
lipid domains (IIlj) . This analysis relies on the above mentioned mutual alignment 
of like domains. Domain-alignment is, however, also of biological relevance, for 
example in the context of, the immune response, where organization of receptor- 
ligand domains occurs during T-cell adhesion mm- Both, the formation of such 
domains as well as the adhesion affinity depend strongly on thermal fluctuations 
and consequently on the bending rigidity of membranes CES1II2D- It is therefore 
reasonable to expect that fundamental intermembrane interactions will play an 
important role also in receptor-ligand domain alignment. 

Within the broad DLVO paradigm (1T51) the fundamental long-range interac¬ 
tions between soft material interfaces, mediated by their molecular environment, 
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such as solvation (hydration) interaction, electrostatic interaction, and van der 
Waals interaction, can be treated independently and additively. However, this 
additivity Ansatz is in general not vindicated for entropically driven bending un¬ 
dulation interactions that warrant a more sophisticated approach (I18H20D . 

Besides the fundamental role of entropic membrane undulations, their relation 
with the membrane bending rigidity K c dHD, and through it their connection with 
diverse physiological processes, has spurred a sustained scientific interest m 
Shape analysis of giant unilamellar vesicles (GUV) (1221) . diffuse X-ray scattering 
from oriented lipid multibilayers (1231) . and GUV micropipette aspiration (1241) are 
all techniques exploiting this connection, but none of them so far has been able 
to simultaneously determine the bending moduli for coexisting membrane phases. 
On the other hand, macroscopically sized domains form distinct lamellar lattices 
in multibilayer systems, enabling the application of osmotic stress experiments 
(PESO). In such experiments, osmotic pressure is maintained by, e.g., large neu¬ 
tral polymers, such as poly-ethylene-glycol (PEG), which do not penetrate into 
the interbilayer water layer, while the corresponding bilayer separation and more 
recently also the specific line broadening due to fluctuations are measured by small- 
angle X-ray scattering (SAXS). Several groups, including ours, have previously ap¬ 
plied this approach to study interactions between macromolecules, including lipid 

bilayers (JHl G5H331) • 

The bare long-range DLVO interaction components, that couple macromolec- 
ular surfaces through their molecular environment, get inextricably intertwined 
through the thermally driven conformational fluctuations of the soft interfaces, 
making detailed predictions of the overall interaction nearly impossible. There- 
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fore, many studies in the past have resorted to describe such complicated ther¬ 
mal fluctuation effects by different mean-field/additivity approximations, where 
conformational fluctuation effects on the bare interaction potentials are included 
self-consistently (jT3} 1201 135H371) . In contrast, additivity/mean-field approximations 
can be altogether avoided in the case of simulations that start from fundamental 
long-range DLVO interaction components and need no additional approximations 
to yield an accurate estimate for the total osmotic pressure in the system (1381 132} • 
In order to understand the coupling between bare interactions and thermal 
undulations in phase separated systems, we apply a gradient-based optimization 
algorithm to iteratively adjust the parameters entering MC simulations, i.e., the 
coefficients describing the strength and range of intermembrane interactions as 
well as the bending rigidity characterizing the thermal undulations, in order to 
best match simulation results with the experimental osmotic stress data for co¬ 
existing Lo/Ld phases. We demonstrate the capability of the simulation-driven 
analysis choosing a well-studied mixture of dioleoyl phosphatidylcholine (DOPC), 
distearoyl phosphatidylcholine (DSPC) and cholesterol (Choi) (130H32)) . previously 
shown to exhibit Lo/Ld domain alignment in the phase coexistence regime (1111) . 
We find that Lo domains are about three times more rigid than Ld domains, which 
exhibit significant contributions to domain repulsion from bending fluctuations. 
On the other hand, hydration forces decay much slower with domain separation 
between Lo domains. In turn, attractive van der Waals interactions were found to 
be of similar magnitude between Lo-Lo and between Ld-Ld domains. Our results 
provide insight into the strength and distance dependence of forces at play between 
like-domains as a prerequisite to devising theories for domain alignment. 
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Materials and methods 

Sample preparation 

DSPC, DOPC, and Choi were purchased from Avanti Polar Lipids, Inc., Alabaster, 
AL, USA and used without further purification. Poly(ethylene glycol) (PEG) with 
an average MW of 8000 was obtained from Fluka Chemie AG, Buchs, Switzerland 
and used as received. 

After weighing, lipids were dissolved in chlorofornr/methanol 2:1 at concen¬ 
trations of 10 mg ml -1 (14311 . We prepared the ternary lipid-only mixture DOPC/ 
DSPC/Chol (0.42:0.37:0.21) in a glass vial and evaporated the organic solvent 
under a gentle nitrogen stream at 30 °C. This lipid composition and its tie-line 
lie well inside the Lo/Ld phase coexistence region according to (14 0[ I4TD . and the 
domains’ structural properties have already been investigated with different meth¬ 
ods dm as}. Remaining solvent traces were removed by placing the samples in 
vacuum overnight. 18 Mil cm -1 water (UHQ PS, USF Elga, Wycombe, UK) was 
added at 20 pi water /mg lipid and the mixtures fully hydrated at 50 °C for 4 hours 
with repeated freeze-thaw cycles. 

To exert osmotic pressure on MLVs, samples were cooled to room temperature 
after hydration and aliquots overlayed with PEG dissolved in water, yielding final 
concentrations of 1-42 wt% PEG in water. Samples were protected against oxida¬ 
tion with argon, the vials closed and taped, and stored at 4 °C for 7-10 days until 
the measurement. The osmotic equation of state for PEG, connecting its osmotic 
pressure with its solution concentration is well known (14411 and allows for an accu¬ 
rate determination of the PEG osmotic pressure P by using previously reported 
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high resolution data (145|) . 

X-ray measurements 

Small-angle X-ray scattering (SAXS) was performed at the Austrian SAXS beam¬ 
line at ELETTRA, Trieste, Italy (1461147jl . at a wavelength of 1.54 A and an energy- 
dispersion A E/E of 2.5 x ICE 3 . A mar300 Image Plate 2D detector from mar- 
research, Norderstedt, Germany was used, covering a g-range from 0.2-7.1 A -1 
and calibrated with silver-behenate (CH 3 (CH 9 ) 20 —COOAg) with a d-spacing of 
5.838 nm (|48l) . Samples were filled into reusable quartz-glass capillaries and kept 
in a brass sample holder connected to a circulating water bath from Huber, Offen- 
burg, Germany. The samples were equilibrated for 10 min at (20.0 ± 0.1) °C before 
exposing them for 30 s to the X-ray beam. 

The two dimensional detector signal was radially integrated with FIT2D (1491 
EDI). Water background subtraction for samples without PEG was performed with 
Primus (I5ip . For osmotically stressed samples however, additional scattering from 
PEG made a standard background subtraction impractical. Since the essential 
informations in this case were just the Bragg peaks’ shapes and positions, we 
subtracted approximative backgrounds, obtained by interpolating between SAXS 
signals of pure water and PEG/water mixtures. Alternatively, one could just 
subtract an arbitrary smooth function from the measured data. 

The reduced data were then fitted using a recently published, full g-range anal¬ 
ysis method for coexisting liquid/liquid membrane domains (1521) . This method 
models each phase’s contribution individually with a bilayer-structure and a su¬ 
perimposed membrane lattice. The lattice description is based on a modified Caille 
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theory (15511531) and therefore yields the average membrane periodicity d and the 
line shape parameter 77 , which is connected to the mean square fluctuation of the 
membrane spacing via A 2 = rjd 2 /ir 2 (132 j) . The bilayer-structure of each phase is 
then modeled separately via probability distributions of quasi-molecular fragments 

(USD- 

Most importantly, the full g-range analysis allowed us to quantify the magni¬ 
tude of fluctuations for coexisting domains. For both phases of stress-free sam¬ 
ples, this also yields accurate electron density profiles, from which the bilayer 
thickness could be obtained; but this was not possible when osmotic pressure was 
applied. Instead, the osmotic thickening of ds was calculated using ds{P) = 
ds( 0) • {Ka + P ■ d{P ))/ (Ka + P ■ ds( 0)) (flTD) . where the area extension modulus 
Ka was estimated from published micropipette aspiration experiments on single 
lipids and binary lipid mixtures (j56l 1571). as detailed in Sec. SI of the Supporting 
Material. The overall analysis was rather insensitive to uncertainties in I\a be¬ 
cause the maximal change in bilayer thickness was only slightly larger than the 
uncertainty of the fit (±2%). The definition of the bilayer thickness ds was found 
to be more important. In principle one could determine optimal values of ds via 
a joint fit with free MC parameters, but this problem is under-determined and led 
to bizarre values of ds for different data sets (j58|) . Instead, we defined ds as the 
distance between the remotest lipid atoms (|59|) . also known as the steric bilayer 
thickness d22D; this yielded good fits and comparable results, while being directly 
accessible from the SAXS analysis. 
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Figure 1: Real space snapshots of equilibrated Ld simulations at zero (left) 
and finite (5.5 MPa, right) osmotic pressure. Membranes are drawn with 
their average thickness. Deviations from the periodic lattice are color coded. 
Due to 3D periodic boundary conditions, the top-most (orange) and bottom¬ 
most membranes are equal. The most prominent effects of external pressure, 
a compression of the stack and a reduction of the fluctuations, are clearly 
visible. 

Membrane Monte-Carlo simulation 

The simulation code used has been described previously in detail for a single 
membrane between two walls and for a stack of membranes (1381 1581) . For 

completeness, but also to highlight our modifications, we briefly summarize its 
basic elements. 

The system under consideration consists of a stack of M fluctuating and in¬ 
teracting membranes of size L x L, as depicted in Fig. [l] The m-th membrane’s 
displacement from its average plane is denoted as u m (x,y), the average distance 
between membranes a, and the bending rigidity as K c . Imposing periodic bound¬ 
ary conditions in ah directions yields the Hamiltonian of a stack of membranes 

M—l / TV" \ 

Li = ^2 I i^(y 2 u m ) 2 + ^(a m )jdxdy, 

m =0 ' ^ ' 


(1) 
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where $ denotes the bare interaction potential, given here by the hydration repul¬ 
sion and the van der Waals attraction, and a m (x, y) = u m +i(x, y) — u m (x, y)+a de¬ 
notes the local distance between two membranes. We furthermore require a m > 0, 
meaning that membranes cannot interpenetrate. 

To reduce the system’s degrees of freedom to a finite amount, the membranes 
are discretized on a square N x N lattice. The simulation is performed in the 
constant pressure ensemble. (I60|) which converges for this model faster than con¬ 
stant volume simulations (139|) . Monte-Carlo updates are proposed in a and in the 
complex coefficients u m (q x ,q y ) of the Fourier transformation of u m (x,y). Simu¬ 
lating in Fourier space allows for larger moves, thereby accelerating equilibration 
022D- After every Monte Carlo step (MCS), which corresponds to degree of free¬ 
dom (. N 2 M +1) update proposals, we re-centered the coordinate system to correct 
for small center of mass movement as a new feature in the calculations. 

Simulations were performed for L = 700 A, several different N (6, 8, 12, 16, 
24 and 32), M = 8, equilibration lengths of 3 x 10 3 MCS, and collection lengths 
of 10 4 MCS, which typically exceeded the autocorrelation time by a factor of 100. 
Simulations were started with step sizes estimated from an approximative theory 
d2DD and then subsequently optimized during equilibration, applying either dynam¬ 
ically optimized Monte Carlo (DOMC), or - as a new feature - the acceptance ratio 
method (ARM) as a backup if DOMC fails (|58l (UTl) . 

Several observables can be determined from converged simulations, but the 
two most important quantities for comparison with SAXS experiments are the 
temporally- and spatially-averaged distance between membranes dw = (a) and 
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the time average of its fluctuations 

A 2 = ((z m+ i(x,y) - z m (x,y) - d B - d w ) 2 ), (2) 

where the long bar denotes spatial averaging over ( m,x,y ), (.) denotes time aver¬ 
aging, and z m (x, y) = u m (x, y) + m ■ (a + d B ) is the m-th membrane’s position in 
real-space. Specifically, dw corresponds to the experimental thickness of the water 
layers separating the lipid bilayers, while A is related to the experimental Caille 
parameter y as detailed above. 

It should be emphasized that our explicit purpose of making contact with the 
X-ray structure factor and the interactions between bilayers, requires much larger 
systems than can be presently envisioned either for all-atom simulations, used to 
obtain electron density profiles, or even for the most coarse grained molecular sim¬ 
ulations (|62l) . We require M bilayers in a stack, each bilayer having a large lateral 
size L. It has been shown in previous work (j381) . that L = 700 A and M = 8 are 
sufficient to obtain accuracies of 1% for dw and A, and that would require about 
130 000 lipids with associated water in typical molecular simulations. Apart from 
simulation size, also the necessary timescales, which scale with the fourth power of 
the undulation wavelength (|63l pp. 77-78), render molecular dynamics simulations 
for that purpose unfeasible. Furthermore, to fit the experimental data requires on 
the order of 100 separate simulations, distributed on multiple optimizations from 
different start points. In the membrane MC simulations we employ, each bilayer is 
reduced to a network consisting of N nodes in each of the two lateral directions and 
each node has only one degree of freedom. Computed observables change signifi¬ 
cantly with NjL (1381 f3TTT) , so simulations were performed for a sequence of values 
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of N E {6, 8,12,16,24, 32} and then the observables were extrapolated towards 
N/L —>• oo. Further details of this finite size convergence are given in Sec. S2. 


Bare interaction potentials 

For uncharged membranes, the potential at bilayer separation a is modeled canon¬ 
ically by (IM|) 

®(“) ~ ' 4Aexp ( - x) " 12^?' (3) 

The first term is the well-established empirical form of the solvent-mediated hy¬ 
dration interaction, which has been argued to originate from changes in various 
measures of order for the water structure at the membrane interface (I65H67D , with 
the strength A and the decay length A, which is typically in the range of 1-2 A 
(13211 . The second term describes the ubiquitous van der Waals interaction poten¬ 
tial for two planar semi-infinite layers, with H being the Hamaker coefficient that 
in general also depends on the bilayer separation a, H = H(a) (|68l p. 15). This 
functional form is convenient because it can in fact describe both cases of either 
two finite-thickness layers interacting across a solvent layer (|69D , as well as effective 
pairwise interactions in an infinite stack of finite-thickness layers m- For large 
solvent layer thickness the nonpairwise additive effects in the latter case become 
negligible and the van der Waals interaction potential for the two cases follows 
exactly the same separation dependence. 

Due to the divergence of the van der Waals potential for a —> 0, the 1/a 2 term is 
cut off for a < 1A (13811 . In experiments, the collapse of charge neutral bilayers due 
to van der Waals forces is avoided by very short range steric interactions established 


by McIntosh et al. (HD, but which occur at significantly higher osmotic pressures 
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Figure 2: Hamaker coefficient H for hydrocarbon multilayers of height ds 
and separation dw in water. Highlighted are the applied values of H for Ld 
and Lo, which are described in the main text. 

than those relevant for the present experiments, see also Fig. S5. 

To calculate the Hamaker coefficient H ab initio , we had to approximate the 
lipid bilayers by pure hydrocarbon CSJ- Specifically, we calculated H for an infinite 
stack of hydrocarbon layers in water, based on a full multilayer Lifshitz formulation 
Q7Uj). The ranges for the hydrocarbon thicknesses ds = 45 A to 60 A and the 
water spacings d\y = 5 A to 30 A were motivated by our experimental data. In 
this calculation range, differences in the Hamaker coefficient were within 10%. 
For our MC simulations the exact value of H matters most when all forces are 
of comparable magnitude, that is at vanishing external osmotic pressure. We 
therefore used the H values of 4.08 x 10~ 21 J = 4.08 zJ for Ld and 4.15 zJ for Lo 
domains, see Fig. [2] 

Both components of the bare potential, i.e. hydration and van der Waals, cause 
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partial bare pressures between neighboring membranes given by C2D 


Phyd(dw) = A exp 



Pvdw{dw) — ~ Q nd 3 • 


(4) 


For comparison to previous reports using mean-field/additivity approximations 
for modeling undulation interactions, one can obtain an effective decay constant 
A und by subtracting bare contributions from experimental data, i.e. P un d = P — 
Phyd — PvdW ([22) • The undulation decay constant then results from a fit of 
Pund = Aund exp(-A und /dw), with the two adjustable parameters A und and \ und . 
Because the undulation pressure deviated significantly from a perfect exponential, 
we limited the fit to large separations (d\y > 14A). 

Optimizing parameters against experimental data 

Calculation of the Hamaker coefficient H, as described above, allowed us to reduce 
the number of free fitting parameters for the simulations to three, A = (A, A, K c ), 
for a joint analysis of domain separation and fluctuation data (see below). 

We implemented a least squares routine with Matlab® (1741) . utilizing its trust 
region reflective optimization algorithm to minimize the sum of the squared residues 



2 


(5) 


where d\y,i and Aj are the experimentally determined values at fixed osmotic pres¬ 
sure P{, dw(Pf, A) and A(Pf, A) are simulation results, and U e g(f) is the effective 
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uncertainty of a given quantity /, derived from 

£4ff(/) = u 2 (f ex p ) + u 2 (f sim ) + • c/(po) 2 ■ (6) 

The agreement between model and data was evaluated by the reduced Xred = 
X 2 /N, where N equals the number of data points minus the number of free pa¬ 
rameters m p. 268). The Jacobian for this gradient based algorithm and the 
derivative in Eq. ^ were computed with the histogram reweighting method de¬ 
scribed in Sec. S3. Once the iteration converged, the uncertainties of the fit pa¬ 
rameters were determined from the curvature of Xred • 01 'der to locate the global 
optimum, several iterations from randomly chosen initial parameter sets were per¬ 
formed. 

To test our implementation, we fitted simulation results determined for one 
reasonable parameter set A', by starting the least squares from several different 
initial starting points A. Within 3-5 iterations, these optimizations converged 
towards the correct values A', thereby indicating that the weighted histogram based 
differentiation and the fit were correctly implemented. For the experimental data 
sets, convergence was usually reached within 10 iterations. However, due to the 
stochastic nature of the simulations and the consequential randomness of results 
and derivatives, the optimization algorithm propagated poorly in flat regions, i.e. 
small V Xred ■ Because X^ erf (A) is a smooth function and its gradient has to vanish 
at extrema, the optimization algorithm’s efficiency decreased, the closer it got to 
the optimum. This was another reason for starting several independent iterations 

(USD- 

As a further test case, we re-analyzed previously published osmotic pressure 
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data of pure dimyristoyl-phosphocholine (DMPC) bilayers (132 p . yielding very rea¬ 
sonable values and a good agreement between simulations and experiments. De¬ 
tails are given in Sec. S4. Thus, we conclude that our method provides a robust 
analysis for interactions in fluctuating membrane assemblies. 


Results and Discussion 

X-ray analysis 

SAXS patterns were analyzed as detailed previously by a Caille theory-based anal¬ 
ysis (111)) . Figure [3] showcases the analysis for two samples at osmotic pressures of 
34kPa and 2.4 MPa, demonstrating that shapes and positions of Bragg reflections 
are well reproduced. Consistent with previous studies (0131 DU, we find sharper 
and more prominent Bragg reflections for the Lo phase due to its decreased bend¬ 
ing fluctuations, compared to the coexisting Ld phase. Fits for all other samples 
are shown in Sec. S5. For increased osmotic pressures, Bragg peaks shifted towards 
higher q and became more prominent. This is due to the decrease of bilayer sep¬ 
aration which goes in hand with a reduction of bending fluctuations in agreement 
with previous reports 4521177) . 

Peak line-shapes for Lo and Ld domains were found to be well described by 
the applied Caille theory, particular at low osmotic pressure (Fig. S4). Since this 
theory is incapable of fitting peaks from lamellar gel phases (178) , we conclude that 
neither peaks assigned to the Lo, nor to the Ld phase can originate from a gel 
phase. This is also consistent with reported compositional DSPC/DOPC/Chol 
phase diagrams mm and a recent SAXS study from our lab, reporting for the 
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Figure 3: Calculated scattering intensities (solid lines) from full q -range anal¬ 
yses, compared with recorded SAXS data from coexisting phases (dots) for 
two different osmotic pressures P. Bragg reflections from aligned Lo and Ld 
domains are indicated by symbols O and X, respectively. 

identical lipid mixture that the structural parameters match those of Lo and Ld 
phases at the tie-line endpoints (HD- 

Fit quality of SAXS spectra worsened for increased PEG concentrations, see 
Fig- m or Sec. S5. Probably the underlying Caille theory loses its applicability 
for the increased order experienced at elevated osmotic pressures. While effects 
on domain separation were negligible, fluctuations determined from the fits be¬ 
came increasingly skewed with osmotic pressure, in particular for Lo domains (see 
below). 

The effect of osmotic pressure on the lamellar repeat spacing d, as determined 
from the SAXS analysis, is plotted in Fig. [4j At high osmotic stress, the distance 
between bilayers is effectively set by the repulsive hydration interaction which 
dominates the repulsive fluctuation interaction and the attractive van der Waals 
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Figure 4: Osmotic pressure P vs membrane periodicity d for Ld and Lo 
determined by SAXS analysis (1731) . Dashed lines are meant solely as a guide 
for the eye. 

interaction. As osmotic pressure is decreased, the water spacing between bilayers 
dw increases and the fluctuation interaction eventually dominates the hydration 
interaction. As the osmotic pressure is reduced to zero, the attractive van der 
Waals force balances the total repulsive forces, resulting in finite dw and d values. 

Within experimental uncertainty, the two isotherms in Fig. [4] are rather sim¬ 
ilar when the difference in membrane thickness is taken into account (dg d = 
(48.5 ±1.0) A and djf = (61.3 ± 1.2) A). Of course, identical isotherms would 
imply that all the interactions are identical. However, significant experimental 
differences were observed in the fluctuation behavior as detailed below, corrobo¬ 
rating the crucial advantage of jointly analyzing fluctuations and osmotic pressure 
isotherms in order to obtain the interaction parameters (l32l) . 
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Tabic 1: Optimal parameters determined for describing the coexisting Lo/Ld 
phases in DOPC/DSPC/Chol (0.42:0.37:0.21). Errors as obtained from the 
fitting routine, see text for further details. 



Ld 

Lo 

K c / z J 

44 ± 10 

120 ± 20 

A/ Pa 

-^ q 8.3 ± 0.2 

1 Q 8.1 ± 0.2 

a/A 

1.37 ±0.15 

1.74 ±0.15 

X.red 

1.5 ±0.5 

5.8 ±0.5 


Optimized simulations 

The experimental data and the results of optimized simulations are compared in 
Fig- HI while Tab.[l]lists results for the interaction parameters. Experimental errors 
for dw and ij were obtained from the SAXS analysis and for P were estimated 
to equal the pipetting error of 6% for viscous PEG solutions. To quantify the 
agreement between data and simulations, we report Xred> which becomes ca. 1 if 
the differences are compatible with experimental errors (1751 p. 268). This is the 
case for the Ld phase, where simulations and experimental data match ideally, but 
the mismatch for Lo is bigger than expected (xfed. = 6)- 

We are inclined to attribute this discrepancy for Lo at least partially to the 
limited applicability of the Caille theory for highly ordered systems, as described in 
the previous section. Indeed, deviations in A are especially pronounced for small 
bilayer separations, i.e. at high osmotic pressures. In light of these discrepancies, 
we suggest that the experimental uncertainties determined for the Lo phase are 
rather too small because they do not take into account the decreasing applicability 
of the Caille theory for more ordered phases whose fluctuations are suppressed by 
low hydration. 

While differences in P(d\y) are insignificant between Ld and Lo (see also 
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d w (A) d w (A) 


Figure 5: Osmotic pressure (top) and fluctuations (bottom) vs water-layer 
thickness for best fit of membrane MC simulation (cyan/orange) against 
SAXS data (gray) ffTiJD . Solid lines were obtained by exponentially interpo¬ 
lating fluctuation contributions. 
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Fig. 0, fluctuations of the Lo phase, containing most of the DSPC and about 
thrice as much cholesterol as Ld, are evidently smaller (Fig. 0. In the continuum 
mechanics treatment used in the simulations, this increase in bilayer stiffness is 
captured by a threefold higher K c for Lo, see Tab. |Tj 

The values obtained by us for K c compare well with previously reported results 
from different techniques. Bending rigidities of binary DOPC/cholesterol mixtures 
have been measured by several groups, ranging from (60 ± 8) to (100 ± 25) zJ and 
were found to be largely unchanged by the cholesterol content (l80fl83p . This sup¬ 
ports the K c = (44 ± 10) zJ obtained for Ld, where DOPC is the main constitutent 
m>- In contrast, a larger concentration of saturated lipids, for which K c does in¬ 
crease with cholesterol <m, is present in the Lo phase, so a larger bending rigidity 
would be expected for Lo than for Ld. Our finding of K c = (120 ± 20) zJ for the 
Lo phase is consistent with this expectation. 

Furthermore, molecular dynamics (MD) simulation results are available for 
comparison. Khelashvili et al. (I84h used the reported tie-line endpoint compo¬ 
sitions m to separately simulate the liquid-disordered and -ordered phases, ob¬ 
taining bending moduli of 80-130 zJ for Ld and 340-440 zJ for Lo. Although these 
values are large compared to our results, both methods find a strong increase of 
K c between Ld and Lo. 

In agreement with Ref. [MI we find that a rather simple model suffices to relate 
bending to area extension moduli for cholesterol-rich samples (j86p . Based on the 
assumption that the main contribution to membrane rigidity comes from the stiff 


cholesterol ring of size S\ Pan et al. used the relationship 5 72 = 12K C /Ka- For our 
samples, with Ka = 430mNm^ 1 and 2100mNnr -1 (see Sec. SI for details), this 
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equation yields S' = 1 1 A and 8 A for Ld and Lo, respectively, in good agreement to 
actual cholesterol ring sizes of about 9 A, giving additional support to our analysis. 

Interdomain forces 

As stated before, the differences between Ld and Lo in the P vs dw data sets are 
small. However, a more thorough investigation of these quantities yields interesting 
insights. Because good fits to these data were obtained, the total pressure P is 
readily dissected into its individual contributions from the fundamental surface 
forces, whose functional dependences are plotted in Fig. [6] 

The thicker Lo bilayer causes an increase in the Hamaker coefficient, but only 
by 3% compared to the Ld phase; this is a minor difference in the van der Waals 
interaction that is hardly noticeable in the P v dw curve in Fig. [6} For small bilayer 
separations, the hydration interactions are of similar magnitude and represent, as 
expected, the dominant contribution to the total interaction potential for both 
phases. Despite these similarities, the fluctuation pressure starts to surpass the 
hydration pressure already at much smaller separations dw for Ld than for Lo. 
This difference implies, in contrast to the ordered phase, that the undulation 
interaction becomes the most important repulsive interaction over a wider range 
of bilayer separations in the case of the disordered phase. Stronger repulsions due 
to fluctuation interactions are of course reasonable because thermal undulations 
were found to be significantly increased for the Ld phase (Fig. [5j. Nevertheless, 
even in the Lo phase, the thermal undulation interaction dominates the hydration 
force over the most important, well hydrated range of dw, starting at separations 
of 12 A. 
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dvj( A) 


Figure 6: Partitioning of total pressure P into contributions from hydration 
Phydi van der Waals P v dw, and undulations P un d for Ld (top) and Lo (bottom) 
(1791) . The large open black circles show the values of the separation dw, at 
which hydration and undulation pressure are equal. Due to the additive 
constant 5, the hydration pressure deviates from a straight line at low P. 
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We obtained almost exponentially decreasing fluctuation forces of the scaling 
form oc exp (— z/X un d), with effective decay lengths \ U nd ~ 3.3 A and 3.7 A for Ld 
and Lo, respectively The ratio of fluctuation to hydration decay length \ un d /A is 
obtained as 2.4 for Ld and 2.1 for Lo. While the mean-field theory predicted its 
value as 2. o m, values of 2.4 have been reported for simulations (1351 1551) . and 
2.0-3.0 from other experiments (j8j [32| 1331) . 

Compared to Lo, a significantly shorter decay length for the hydration interac¬ 
tion pressure was found for the Ld phase. At present, the origin for this difference 
is unclear. However, it is this difference combined with the larger fluctuation force 
that gives P versus d\y curves that are nearly the same for Lo and Ld, both with 
fully hydrated dw close to 17 A. 

Domain alignment across interlamellar aqueous phases has recently been hy¬ 
pothesized to be caused by water network mismatch due to the different hydration 
properties of Lo and Ld phases ([3]). In support of this postulation, we observed 
significantly different hydration forces and nearly equal van der Waals forces for 
both phases. Thermal fluctuations were however neglected in the aforementioned 
hypothesis, while we now find considerable differences specifically in the undula¬ 
tion forces for coexisting domains. Their importance is especially striking near 
full hydration, where undulation and van der Waals pressures surpass hydration 
repulsion by an order of magnitude (see Fig. [6]). 


Conclusion 

We have evaluated the fundamental long-range interactions between coexisting Lo 
and in Ld domains in DOPC/DSPC/cholesterol, which is a frequently used model 
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system for mammalian outer plasma membranes flnHEg SMI E2 ESD- Because 
we could do this at concentrations where Lo and Ld domains coexist, we were able 
to avoid all uncertainties in the phase diagram and its associated tie-lines between 
Lo and Ld phases. This work combines methodology from three separate inputs: 
SAXS/osmotic stress experiments, comprehensive Monte Carlo simulations, and 
detailed calculations of van der Waals interactions. 

The reported values for fundamental surface forces and bending moduli are the 
first of their kind being, directly obtained from coexisting Lo/Ld domains. The 
underlying full g-range SAXS analysis allowed us to quantify the extent of fluctua¬ 
tions and capture their dependence on osmotic pressure, which proved essential for 
determining the bending rigidities of cholesterol-rich phases. We obtained bend¬ 
ing moduli of 44 zJ for Ld and a roughly threefold higher value for Lo domains, 
attributable to their larger concentrations of saturated lipid and cholesterol. 

While we obtained almost identical van der Waals interactions for aligned Lo 
and Ld domains, the remaining interactions, however, turned out to be strikingly 
different: decay lengths of the hydration pressures differed by 25% between Lo 
and Ld phases, and repulsions due to thermal fluctuations were found to be sig¬ 
nificantly increased for Ld. These findings provide evidence that a combination of 
hydration repulsion and the fluctuation-driven undulation repulsion must be con¬ 
sidered in any quantitative explanation of the long-range positional correlations 
between aligned Lo and Ld domains. In particular the strong entropic contribu¬ 
tion from undulating Ld domains may be a leading term to be considered. We 
therefore expect that our study will form the base for a concise theory of domain 
alignment. 


Lipid domain interactions 


25 


Author Contributions 

B.K. designed and performed research, analyzed data and wrote the paper; P.H. 
designed and performed research and analyzed data; R.P. and J.F.N. contributed 
analytic tools and wrote the paper; G.P. designed and performed research and 
wrote the paper. 


Acknowledgments 

This work is supported by the Austrian Science Fund FWF, Project no. P24459- 
B20 to GP. The computational results presented have been achieved using the 
Vienna Scientific Cluster (VSC). The authors thank Alexander Rieder and Heinz 
Amenitsch for experimental support and Hans Gerd Evertz for critical review of the 
simulation and advice regarding finite size convergence. Support for the original 
development of the MC software was provided to JFN under grant GM44976 from 
the U.S. National Institutes of Health. RP would like to acknowledge the SLO-A 
bilateral grant Nl-0019 of the Slovene Research Agency. 


Supporting Citations 


References (1MWIU3I) appear in the Supporting Material. 



Lipid domain interactions 


26 


References 

1. Semrau, S., T. Idema, T. Schmidt, and C. Storm, 2009. Membrane-Mediated 
Interactions Measured Using Membrane Domains. Biophys. J. 96:4906-4915. 

2. Ursell, T. S., W. S. Klug, and R. Phillips, 2009. Morphology and interaction 
between lipid domains. PNAS 106:13301-13306. 

3. Tayebi, L., Y. Ma, D. Vashaee, G. Chen, S. K. Sinha, and A. N. Parikh, 
2012. Long-range interlayer alignment of intralayer domains in stacked lipid 
bilayers. Nat. Mater. 11:1074-1080. 

4. Karmakar, S., and V. A. Raghunathan, 2005. Structure of phospholipid- 
cholesterol membranes: An x-ray diffraction study. Phys. Rev. E 71:061924. 

5. Chen, L., Z. Yu, and P. J. Quinn, 2007. The partition of cholesterol be¬ 
tween ordered and fluid bilayers of phosphatidylcholine: A synchrotron X-ray 
diffraction study. Biochim. Biophys. Acta, Biomembr. 1768:2873-2881. 

6. Mills, T. T., S. Tristram-Nagle, F. A. Heberle, N. F. Morales, J. Zhao, J. Wu, 
G. E. S. Toombes, J. F. Nagle, and G. W. Feigenson, 2008. Liquid-Liquid 
Domains in Bilayers Detected by Wide Angle X-Ray Scattering. Biophys. J. 
95:682-690. 

7. Staneva, G., C. Chachaty, C. Wolf, K. Koumanov, and P. J. Quinn, 2008. 
The role of sphingomyelin in regulating phase coexistence in complex lipid 
model membranes: Competition between ceramide and cholesterol. Biochim. 
Biophys. Acta, Biomembr. 1778:2727-2739. 



Lipid domain interactions 


27 


8. Pabst, G., B. Boulgaropoulos, E. Gander, B. R. Sarangi, H. Amenitsch, V. A. 
Raghunathan, and P. Laggner, 2009. Effect of Ceramide on Nonraft Proteins. 
J. Membr. Biol. 231:125-132. 

9. Yuan, J., A. Kiss, Y. H. Pramudya, L. T. Nguyen, and L. S. Hirst, 2009. So¬ 
lution synchrotron x-ray diffraction reveals structural details of lipid domains 
in ternary mixtures. Phys. Rev. E 79:031924. 

10. Uppanroochikkal, P., S. Tristram-Nagle, and J. F. Nagle, 2010. Orienta¬ 
tion of Tie-Lines in the Phase Diagram of DOPC/DPPC/Cholesterol Model 
Bionrembranes. Langmuir 26:17363-17368. 

11. Heftberger, P., B. Kollnritzer, A. A. Rieder, H. Amenitsch, and G. Pabst, 
2015. In Situ Determination of Structure and Fluctuations of Coexisting 
Fluid Membrane Domains. Biophys. J. 108:854-862. 

12. Heberle, F. A., and G. W. Feigenson, 2011. Phase Separation in Lipid Mem¬ 
branes. Cold Spring Harb. Perspect. Biol. 3:a004630. 

13. Marsh, D., 2009. Cholesterol-induced fluid membrane domains: A com¬ 
pendium of lipid-raft ternary phase diagrams. Biochim. Biophys. Acta, 
Biomembr. 1788:2114-2123. 

14. Monks, C. R. F., B. A. Freiberg, H. Kupfer, N. Sciaky, and A. Kupfer, 
1998. Three-dimensional segregation of supranrolecular activation clusters in 
T cells. Nature 395:82-86. 


15. Grakoui, A., S. K. Bromley, C. Sumen, M. M. Davis, A. S. Shaw, P. M. Allen, 



Lipid domain interactions 


and M. L. Dustin, 1999. The Immunological Synapse: A Molecular Machine 
Controlling T Cell Activation. Science 285:221-227. 

16. Rozycki, B., R. Lipowsky, and T. R. Weikl, 2010. Segregation of recep¬ 
tor-ligand complexes in cell adhesion zones: phase diagrams and the role of 
thermal membrane roughness. New J. Phys. 12:095003. 

17. Hu, J., R. Lipowsky, and T. R. Weikl, 2013. Binding constants of membrane- 
anchored receptors and ligands depend strongly on the nanoscale roughness 
of membranes. PNAS 110:15283-15288. 

18. Israelachvili, J. N., 2011. Interactions of Biological Membranes and Struc¬ 
tures. In Intermolecular and surface forces, Academic Press, Burlington, MA, 
577-616. 

19. Helfrich, W., 1978. Steric Interaction of Fluid Membranes in Multilayer 
Systems. Z. Naturforsch., A: Phys. Sci. 33:305. 

20. Podgornik, R., and V. A. Parsegian, 1992. Thermal-mechanical fluctuations 
of fluid membranes in confined geometries: the case of soft confinement. 
Langmuir 8:557-562. 

21. Pabst, G., 2013. Coupling Membrane Elasticity and Structure to Protein 
Function. In Advances in Planar Lipid Bilayers and Liposomes, Elsevier, 
volume 18, 81-109. 

22. Meleard, P., C. Gerbeaud, T. Pott, L. Fernandez-Puente, I. Bivas, M. D. Mi- 
tov, J. Dufourcq, and P. Bothorel, 1997. Bending elasticities of model mem- 



Lipid domain interactions 


29 


branes: influences of temperature and sterol content. Biophys. J. 72:2616- 
2629. 

23. Lyatskaya, Y., Y. Liu, S. Tristram-Nagle, J. Katsaras, and J. F. Nagle, 2000. 
Method for obtaining structure and interactions from oriented lipid bilayers. 
Phys. Rev. E 63:011907. 

24. Evans, E., and W. Rawicz, 1990. Entropy-driven tension and bending elas¬ 
ticity in condensed-fluid membranes. Phys. Rev. Lett. 64:2094-2097. 

25. Boulgaropoulos, B., M. Rappolt, B. Sartori, H. Amenitsch, and G. Pabst, 
2012. Lipid Sorting by Ceramide and the Consequences for Membrane Pro¬ 
teins. Biophys. J. 102:2031-2038. 

26. LeNeveu, D. M., and R. P. Rand, 1977. Measurement and modification of 
forces between lecithin bilayers. Biophys. J. 18:209-230. 

27. Parsegian, V. A., N. Fuller, and R. P. Rand, 1979. Measured work of defor¬ 
mation and repulsion of lecithin bilayers. PNAS 76:2750-2754. 

28. Parsegian, V. A., R. P. Rand, N. L. Fuller, and D. C. Ran, 1986. Osmotic 
stress for the direct measurement of intermolecular forces. In Lester Packer, 
editor, Methods in Enzymology, Academic Press, volume 127 of Biomem¬ 
branes Part O: Protons and Water: Structure and Translocation , 400-416. 

29. McIntosh, T. J., and S. A. Simon, 1986. Hydration force and bilayer defor¬ 
mation: a reevaluation. Biochemistry 25:4058-4066. 



Lipid domain interactions 


30 


30. McIntosh, T. J., and S. A. Simon, 1993. Contributions of hydration and steric 
(entropic) pressures to the interactions between phosphatidylcholine bilayers: 
Experiments with the subgel phase. Biochemistry 32:8374-8384. 

31. Rand, R. P., and V. A. Parsegian, 1989. Hydration Forces Between Phos¬ 
pholipid Bilayers. Biochim. Biophys. Acta 988:351-376. 

32. Petrache, H. I., N. Gouliaev, S. Tristram-Nagle, R. Zhang, R. M. Suter, 
and J. F. Nagle, 1998. Interbilayer interactions from high-resolution x-ray 
scattering. Phys. Rev. E 57:7014-7024. 

33. Pabst, G., S. Danner, R. Podgornik, and J. Katsaras, 2007. Entropy-Driven 
Softening of Fluid Lipid Bilayers by Alamethicin. Langmuir 23:11705-11711. 

34. Pabst, G., N. Kucerka, M.-P. Nieh, M. C. Rheinstadter, and J. Katsaras, 
2010. Applications of neutron and X-ray scattering to the study of biologically 
relevant model membranes. Chem. Phys. Lipids 163:460-479. 

35. Sornette, D., and N. Ostrowsky, 1986. Importance of membrane fluidity on 
bilayer interactions. J. Chem. Phys. 84:4062-4067. 

36. Evans, E. A., and V. A. Parsegian, 1986. Thermal-mechanical fluctuations 
enhance repulsion between bimolecular layers. PNAS 83:7132-7136. 

37. Mecke, K. R., T. Charitat, and F. Graner, 2003. Fluctuating Lipid Bilayer in 
an Arbitrary Potential: Theory and Experimental Determination of Bending 
Rigidity. Langmuir 19:2080-2087. 



Lipid domain interactions 


31 


38. Gouliaev, N., and J. F. Nagle, 1998. Simulations of Interacting Membranes 
in the Soft Confinement Regime. Phys. Rev. Lett. 81:2610-2613. 

39. Gouliaev, N., and J. F. Nagle, 1998. Simulations of a single membrane be¬ 
tween two walls using a Monte Carlo method. Phys. Rev. E 58:881-888. 

40. Zhao, J., J. Wu, F. A. Heberle, T. T. Mills, P. Klawitter, G. Huang, 
G. Costanza, and G. W. Feigenson, 2007. Phase studies of model biomern- 
branes: Complex behavior of DSPC/DOPC/Cholesterol. Biochim. Biophys. 
Acta, Biomembr. 1768:2764-2776. 

41. Heberle, F. A., J. Wu, S. L. Goh, R. S. Petruzielo, and G. W. Feigenson, 
2010. Comparison of Three Ternary Lipid Bilayer Mixtures: FRET and ESR 
Reveal Nanodomains. Biophys. J. 99:3309-3318. 

42. Heberle, F. A., R. S. Petruzielo, J. Pan, P. Drazba, N. Kucerka, R. F. Stan- 
daert, G. W. Feigenson, and J. Katsaras, 2013. Bilayer Thickness Mismatch 
Controls Domain Size in Model Membranes. J. Am. Chem. Soc. 135:6853- 
6859. 

43. Supplier provided M\y, accounting for an additional water molecule with 
DOPC, were used for determining stock concentrations. 

44. Cohen, J. A., R. Podgornik, P. L. Hansen, and V. A. Parsegian, 2009. A 
Phenomenological One-Parameter Equation of State for Osmotic Pressures 
of PEG and Other Neutral Flexible Polymers in Good Solventsf. J. Phys. 


Chem. B 113:3709-3714. 



Lipid domain interactions 


32 


45. Stanley, C. B., and H. H. Strey, 2003. Measuring Osmotic Pressure of 
Poly(ethylene glycol) Solutions by Sedimentation Equilibrium Ultracentrifu¬ 
gation. Macromolecules 36:6888-6893. 

46. Amenitsch, H., M. Rappolt, M. Kriechbaum, H. Mio, P. Laggner, and 
S. Bernstorff, 1998. First performance assessment of the small-angle X-ray 
scattering beamline at ELETTRA. J. Synchrotron Radiat. 5:506-508. 

47. Bernstorff, S., H. Amenitsch, and P. Laggner, 1998. High-Throughput Asym¬ 
metric Double-Crystal Monochromator of the SAXS Beamline at ELETTRA. 
J. Synchrotron Radiat. 5:1215-1221. 

48. Huang, T. C., H. Toraya, T. N. Blanton, and Y. Wu, 1993. X-ray powder 
diffraction analysis of silver behenate, a possible low-angle diffraction stan¬ 
dard. J. Appl. Crystallogr. 26:180-184. 

49. Hammersley, A. P., 1997. FIT2D: an introduction and overview. European 
Synchrotron Radiation Facility Internal Report ESRF97HA02T . 

50. Hammersley, A. P., S. O. Svensson, M. Hanfland, A. N. Fitch, and D. Hauser- 
mann, 1996. Two-dimensional detector software: From real detector to ide¬ 
alised image or two-theta scan. High Pressure Res. 14:235-248. 

51. Konarev, P. V., V. V. Volkov, A. V. Sokolova, M. H. J. Koch, and D. I. Sver- 
gun, 2003. PRIMUS: a Windows PC-based system for small-angle scattering 
data analysis. J. Appl. Crystallogr. 36:1277-1282. 

52. Ref. [TTJ We checked the X-ray analysis for coexisting phases by comparing 
it with PEG-free, homogeneous samples prepared at the published tie-line 


Lipid domain interactions 


33 


endpoint concentrations of (0.79:0.09:0.12) for Ld and (0.05:0.65:0.30) for the 
Lo phase <@D- These samples were also helpful for constraining some model 
details (the widths and distances between molecular subgroups composing 
the lipid heads) in the X-ray analysis. 

53. Zhang, R., R. M. Suter, and J. F. Nagle, 1994. Theory of the structure factor 
of lipid bilayers. Phys. Rev. E 50:5047-5060. 

54. Caille, A., 1972. Physique cristalline: remarques sur la diffusion des rayons 
X dans les smectiques. C. R. Acad. Sc. Paris Sie B. 274:891-893. 

55. Kucerka, N., J. F. Nagle, J. N. Sachs, S. E. Feller, J. Pencer, A. Jackson, and 
J. Katsaras, 2008. Lipid Bilayer Structure Determined by the Simultaneous 
Analysis of Neutron and X-Ray Scattering Data. Biophys. J. 95:2356-2367. 

56. Rawicz, W., B. Smith, T. McIntosh, S. Simon, and E. Evans, 2008. Elas¬ 
ticity, Strength, and Water Permeability of Bilayers that Contain Raft 
Microdomain-Forming Lipids. Biophys. J. 94:4725-4736. 

57. Rawicz, W., K. C. Olbrich, T. McIntosh, D. Needham, and E. Evans, 2000. 
Effect of Chain Length and Unsaturation on Elasticity of Lipid Bilayers. 
Biophys. J. 79:328-339. 

58. Gouliaev, N., 1998. Monte-Carlo simulations of membrane systems. Ph.D. 
thesis, Carnegie Mellon University, Pittsburgh, Pennsylvania, USA. 

59. Specifically, we use d B = 2{z C hoiCH3 + °Choicm ), where z C hoiCH3 and 
&ChoiCH3 are the position (measured from the bilayer center) and the width, 


Lipid domain interactions 


34 


respectively, of the lipid head choline’s CH 3 groups. Within measurement 
accuracy, the definition used in Ref. 1521 yields equal values. 

60. McDonald, I., 1972. NpT-ensemble Monte Carlo calculations for binary liquid 
mixtures. Mol. Phys. 23:41-58. 

61. Bouzida, D., S. Kumar, and R. H. Swendsen, 1992. Efficient Monte Carlo 
methods for the computer simulation of biological molecules. Phys. Rev. A 
45:8894-8901. 

62. Cooke, I. R., and M. Deserno, 2005. Solvent-free model for self-assembling 
fluid bilayer membranes: Stabilization of the fluid phase based on broad 
attractive tail potentials. J. Chem. Phys. 123:224710. 

63. Pabst, G., N. Kucerka, M.-P. Nieh, and J. Katsaras, 2014. Liposomes, Lipid 
Bilayers and Model Membranes: From Basic Research to Application. CRC 
Press. 

64. An additional steric repulsion A s t\ s t exp(— a/X s t), with A s t = 3.6 GPa and 
A s t = 0.6 A according to (171)1 was used, but proved unimportant for realistic 
parameters. 

65. Marcelja, S., and N. Radic, 1976. Repulsion of interfaces due to boundary 
water. Chem. Phys. Lett. 42:129-130. 

66. Kanduc, M., E. Schneck, and R. R. Netz, 2013. Hydration interaction be¬ 
tween phospholipid membranes: insight into different measurement ensembles 
from atomistic molecular dynamics simulations. Langmuir 29:9126-9137. 


Lipid domain interactions 


35 


67. Kanduc, M., A. Schlaich, E. Schneck, and R. R. Netz, 2014. Hydration repul¬ 
sion between membranes and polar surfaces: Simulation approaches versus 
continuum theories. Adv. Colloid Interface Sci. 208:142-152. 

68. Parsegian, V. A., 2006. Van der Waals forces. Cambridge university press 
Cambridge etc. 

69. The van der Waals potential for infinitely extended slabs of finite thickness 
ds at separation a is oc (1/a 2 — 2/(a + ds) 2 + l/(a + 2de) 2 ) (1551 p. 15), and 
was applied for osmotic stress experiments in e.g. Refs ® [32, 

70. Podgornik, R., R. H. French, and V. A. Parsegian, 2006. Nonadditivity in 
van der Waals interactions within multilayers. J. Chem. Phys. 124:044709. 

71. McIntosh, T. J., A. D. Magid, and S. A. Simon, 1987. Steric repulsion 
between phosphatidylcholine bilayers. Biochemistry 26:7325-7332. 

72. Although this model gives only a first order estimate for the van der Waals 
interactions of fluctuating lipid bilayers, it is to our knowledge the best avail¬ 
able approximation in the absence of data on the dielectric response of PC 
lipids. Further effects of, e.g. lipid headgroup dipolar moment fluctuations 
(1551) . could be considered as well, but they would be only important at very 
small separations where hydration forces dominate and the exact form of the 
van der Waals interaction is irrelevant. 

73. Equation Q was derived from Pj(dw) ~ — d^j(dw)/ddw- The difference to 
the exact relationship Pj(d\y) = (—d$(a)/da) was found to be less than the 
simulational uncertainty. 


Lipid domain interactions 


36 


74. 2011. MATLAB v. 7.12 (R2011a). 

75. Taylor, J., 1997. Introduction to Error Analysis, the Study of Uncertainties 
in Physical Measurements, 2nd Edition, volume 1. University Science Books, 
New York. 

76. Alternatively, one could have used optimization algorithms specialized for 
simulations (j90H93]), but the existing implementations did not satisfy our 
needs. 

77. Hemmerle, A., L. Malaquin, T. Charitat, S. Lecuyer, G. Fragneto, and 
J. Daillant, 2012. Controlling interactions in supported bilayers from weak 
electrostatic repulsion to high osmotic pressure. PNAS 109:19938-19942. 

78. Pabst, G., 2006. Global properties of biomimetic membranes: perspectives 
on molecular features. Biophys. Rev. Lett. 01:57-84. 

79. The small offset 5 = 200 Pa is necessary for plotting P = 0 on a logarithmic 
scale. 

80. Pan, J., T. T. Mills, S. Tristram-Nagle, and J. F. Nagle, 2008. Cholesterol 
Perturbs Lipid Bilayers Nonuniversally. Phys. Rev. Lett. 100:198103. 

81. Sorre, B., A. Callan-Jones, J.-B. Manneville, P. Nassoy, J.-F. Joanny, 
J. Prost, B. Goud, and P. Bassereau, 2009. Curvature-driven lipid sort¬ 
ing needs proximity to a demixing point and is aided by proteins. PNAS 


106:5622-5626. 


Lipid domain interactions 


37 


82. Tian, A., B. R. Capraro, C. Esposito, and T. Baumgart, 2009. Bending Stiff¬ 
ness Depends on Curvature of Ternary Lipid Mixture Tubular Membranes. 
Biophys. J. 97:1636-1646. 

83. Gracia, R. S., N. Bezlyepkina, R. L. Knorr, R. Lipowsky, and R. Dirnova, 
2010. Effect of cholesterol on the rigidity of saturated and unsaturated mem¬ 
branes: fluctuation and electrodeformation analysis of giant vesicles. Soft 
Matter 6:1472-1482. 

84. Khelashvili, G., B. Kollmitzer, P. Heftberger, G. Pabst, and D. Harries, 2013. 
Calculating the Bending Modulus for Multicomponent Lipid Membranes in 
Different Thermodynamic Phases. J. Chem. Theory Comput. 9:3866-3871. 

85. Pan, J., S. Tristram-Nagle, and J. F. Nagle, 2009. Effect of cholesterol on 
structural and mechanical properties of membranes depends on lipid chain 
saturation. Phys. Rev. E 80:021931. 

86. This relation was suggested by Evan Evans (|851) . 

87. Veatch, S. L., and S. L. Keller, 2003. Separation of Liquid Phases in Giant 
Vesicles of Ternary Mixtures of Phospholipids and Cholesterol. Biophys. J. 
85:3074-3083. 

88. Scherfeld, D., N. Kahya, and P. Schwille, 2003. Lipid Dynamics and Do¬ 
main Formation in Model Membranes Composed of Ternary Mixtures of Un¬ 
saturated and Saturated Phosphatidylcholines and Cholesterol. Biophys. J. 


85:3758-3768. 


Lipid domain interactions 


38 


89. Podgornik, R., 1988. Solvent structure effects in dipole correlation forces. 
Chem. Phys. Lett. 144:503-508. 

90. Ben-Tal, A., and A. Nemirovski, 1999. Robust solutions of uncertain linear 
programs. Operations Research Letters 25:1-13. 

91. Ben-Tal, A., and A. Nemirovski, 2000. Robust solutions of Linear Program¬ 
ming problems contaminated with uncertain data. Math. Program. 88:411— 
424. 

92. Fu, M. C., 2002. Feature Article: Optimization for simulation: Theory vs. 
Practice. INFORMS Journal on Computing 14:192-215. 

93. Ben-Tal, A., A. Goryashko, E. Guslitzer, and A. Nemirovski, 2004. Ad¬ 
justable robust solutions of uncertain linear programs. Math. Program., Ser. 
A 99:351-376. 

94. Salsburg, Z. W., J. D. Jacobson, W. Fickett, and W. W. Wood, 1959. Ap¬ 
plication of the Monte Carlo Method to the Lattice-Gas Model. I. Two- 
Dimensional Triangular Lattice. J. Chem. Phys. 30:65-72. 

95. Ferrenberg, A. M., and R. H. Swendsen, 1988. New Monte Carlo technique 
for studying phase transitions. Phys. Rev. Lett. 61:2635-2638. 

96. Ferrenberg, A. M., and R. H. Swendsen, 1989. Optimized Monte Carlo data 
analysis. Phys. Rev. Lett. 63:1195-1198. 


97. Kumar, S., J. M. Rosenberg, D. Bouzida, R. H. Swendsen, and P. A. Kollman, 



Lipid domain interactions 


39 


1992. The weighted histogram analysis method for free-energy calculations 
on bionrolecules. I. The method. J. Comput. Chem. 13:1011-1021. 

98. Nagle, J. F., 2013. Introductory Lecture: Basic quantities in model biomern- 
branes. Faraday Discuss. 161:11. 

99. Chu, N., N. Kucerka, Y. Liu, S. Tristram-Nagle, and J. F. Nagle, 2005. 
Anomalous swelling of lipid bilayer stacks is caused by softening of the bend¬ 
ing modulus. Phys. Rev. E 71:041904. 

100. Bonner, J. C., and M. E. Fisher, 1964. Linear Magnetic Chains with 
Anisotropic Coupling. Phys. Rev. 135:A640-A658. 

101. Nagle, J. F., and J. C. Bonner, 1970. Numerical studies of the Ising chain 
with long-range ferromagnetic interactions. J. Phys. C: Solid State Phys. 
3:352. 

102. Quenouille, M. H., 1956. Notes on Bias in Estimation. Biometrika 43:353- 
360. 

103. Tukey, J. W., 1958. Bias and confidence in not quite large samples (abstract). 
Ann. Math. Statist. 29:614. 



Supporting information for: 

Bending rigidities and interdomain forces in 
membranes with coexisting lipid domains 

B. Kollmitzer, P. Heftberger, R. Podgornik, 

J. F. Nagle, and G. Pabst 

April 15, 2015, Graz 


SI Area extension modulus estimation 


The dependence of bilayer thickness on osmotic pressure P is accounted for via the area 
extension modulus K A and given by the equation 1 


d B {P) — d B (0) 


I< A + Pd(P) 
K A + Pd B (0)‘ 


(SI) 


We estimated K A for our coexisting liquid phases based on published data for single 
lipids and binary lipid mixtures by Rawicz et al . 2,3 The Ld phase under investigation 
consists essentially of DOPC, with approximately 10mol% cholesterol. 4 Interpolating 
linearly between the two published values for 0 and 50mol% cholesterol in DOPC 2 
yields /Ci(Ld) = (430 ± 30) mN m _1 . 

In the coexisting Lo phase, the main constituent is the saturated lipid DSPC, which 
is accompanied by ca. 30mol% cholesterol. 4 Unfortunately, published Ka values for 
saturated lipids are sparse. As a compromise, we interpolated linearly between pure 
DMPC (0mol% cholesterol) and a 1:1 mixture of sphingomyelin/cholesterol, 2 ’ 3 yielding 
Ka — (2100 ± 500) mN m _1 for our Lo phase. 

As pointed out in the section X-ray measurements of the main text, knowing the 
magnitude of K A is more important than getting the precise number. That is because 
the biggest estimated change in bilayer thickness turned out to be just 0.3 A. In principle, 
such a subtle difference in d B would be resolvable with SAXS, but not with the additional 
scattering signal due to PEG. 


S2 Finite size convergence 

With open edges, one generally expects a ‘surface’ perturbation proportional to the rel¬ 
ative size of the boundary to the interior, i.e. proportional to 1/N for our systems. As is 


SI 



well known, periodic boundary conditions generally reduce this perturbation. They also 
speed up the convergence with system size, from 1 /N to 1 /N 2 in a case well documented 
by Bonner and Fisher 5 (note their Fig. 1) and in the case of the one-dimensional Ising 
model the convergence is exponentially fast with periodic boundary condions. While 
another case with very slow convergence is known, 6 that one is due to very long range 
interactions not present in our membrane stacks. For periodic boundary conditions, the 
exact solution of a harmonic approximation to Eq. (3) suggests that d\\r and A converge 
asymptotically like y(N) ~ — C2/-/V 2 , i.e. convergence is expected to be faster than 

1/N and, in agreement with the previous simulations,' our results are consistent with a 
dominant 1/N 2 asymptotic convergence, allowing, of course, for higher order terms. 

We perform simulations for several ‘densities’ N € { N m i n ,..., N max } and fit them with 
the function y(N) — + J2k=T °k/N k . Together with the originally proposed k max — 

3 and N € {6,...,32}, this method yields sufficiently precise continuum estimators 
Coo, compared to the experimental uncertainties. 8 However, we found that varying the 
arbitrary parameters k max and N m i n influenced the final estimator stronger for some 
simulations (e.g. high pressures) than for others. To obtain more reliable uncertainties 
and perhaps even better continuum estimates, we perform now several extrapolations, 
with different values for k max and N m i n , but always using the highest possible N max . 
By not changing N max , we weight the most significant simulations (with the highest 
density) stronger. This procedure yields a list of results for Coo,j, which we average for 
the final estimator. Its uncertainty is then determined by the individual errors of c^y 
(statistical uncertainty of observables due to finite simulation length) and their standard 
deviation (error due to finite simulation density). This procedure is closely related to 
the Jackknife technique. 9,10 

Comparisons between these improved Jackknife estimators and estimators obtained by 
the original method are given in Fig. SI. The relative difference in the estimators were 
less than 5% for all performed simulations, but most importantly, Jackknife produces a 
meaningful uncertainty. 


S3 Efficient differentiation 


A single simulation of a particular set of parameters A = (P, A, H , A, K c ,...) contains 
more information in the generated time series, than the aforementioned observables 
which are determined by averaging. By reweighting the simulated histogram of density 
of states, it is possible to compute these quantities over a certain range of simulation 
parameters and thereby also derive their gradients. 11-1314 This well recognized method 
was briefly mentioned for membrane MC simulations, 15 but has not been implemented 
for them previously. 

We calculated the expectation value of an observable f(u, d) for a different set of 
parameters A' from a simulation performed at A by 


(f)X' 


E /a/(W a) ■ exp (■ -SG/kT) 
J2 exp (—SG/kT) 


(S2) 
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Figure SI: Finite size convergence of membrane spacing dw vs membrane “density” N 
of Ld domains according to Tab. 1 at intermediate (top) and high osmotic 
pressures P (bottom). A variant of Jackknife allows us to obtain reasonable 
errors for the estimator. Statistical uncertainties for plotted finite N data are 
less than 10~ 2 A. 
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where the sums extend over all realized configurations and SG is the change in the Gibbs 
energy of each state (u, a) upon changing A to A'. Most parameters could be separated 
from u and a in our case, yielding SG(u, a) — SA ■ £(it, a). This allowed us to store only 
the time series of £ instead of all realized states. The parameters P, A, H , and K c were 
separable in this way, yielding 

^ - SP Zp + 6AXU- ^_Zh + N 2SI ^ £ Kc , (S3) 

where = a/d w , £4 = exp(-a/A), — 1/a 2 , and £ Kc — q 4 \u m (q x ,q y )\ 2 . The local 
distance between membranes is denoted by a — u m+ \ ( x , y) — u m (x , y) + a, while the 
bars denote averages over ( m,x,y ) or ( m,q x ,q y ). V — L 2 Ma is the membrane stack’s 
volume. 

Separating A from (it, a) in SG turned out to be impossible, but we were able to 
calculate gradients of dw and A with respect to A efficiently. Because d\v and A did 
not depend explicitly on A (i.e. df/d A = 0), differentiating Eq. (S2) yielded, 

" ~ kTN 2 n ’ ^ S4) 

where sums extend over all realized states, denotes the collection length and 



Up to first order, (f)y was then determined from (f)y ~ (f}\ + (A' — A) d(f)/d A. 

Thus, for any observable / € {dw, A} and parameter A € {P, A, H, K c , A}, we first 
determined (/) i, 2 {N) f° r Ai^ = A±(5A as detailed above, extrapolated these expectation 
values for A" —> oo according to section S2, and finally calculated the finite difference 
quotient d(f)/d A w ((/)i — (f)?)/‘2dK.. Relative finite differences were set to hA/A = 
0.03. 

We checked this method against direct numerical differentiation for a couple of rea¬ 
sonable parameters. Errors were always sufficiently small (well below 50%) to lead the 
optimization routine towards a global minimum (see the section Optimizing parameters 
against experimental data of the main text). 

S4 Results for a homogeneous control sample 

We tested our analysis on already published SAXS data for homogeneous DMPC MLVs 
determined at 30°C. 16 The Lifshitz calculation of the van der Waals forces yielded a 
value of H — 4.11 zJ for the published bilayer thickness of 44.0 A. The obtained values 
describing the intersurface forces are given in Tab. SI, while Fig. S2 compares the sim¬ 
ulations with the experimental data. Reassuringly, the simulations fit the experimental 
osmotic pressure data well. While the fit to A is excellent for high hydration, the fit 
becomes relatively poor for A as dw becomes small, similarly to our L a sample and 
likely for the same reason given in the main text. 


d(f)x> 

OX 1 A'=A 


S4 




4 6 8 10 12 14 16 18 20 

ata (A) 

Figure S2: Osmotic pressure (top) and fluctuations (bottom) vs water-layer thickness 
for best fit of membrane MC simulation (cyan) against SAXS data (light 
gray) obtained from Ref. 16. 17 Solid lines were obtained by exponentially 
interpolating fluctuation contributions. 
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The interaction parameters obtained from the fit are shown in Tab. SI. Literature 
values for DMPC’s bending modulus range from 50-130 zJ at 30°C. 18 In light of this 
large variation, comparing only results of related methods is appropriate. Ref. 16 could 
not determine K c and the modulus B separately and therefore considered several values 
of K c \ two of these are shown in Tab. SI. The values of A agree very well with ours. 
The larger values of A would have been smaller if the true value of Ka had been known 
at that time. Two differences from the previous analyses are that here we calculated 
H and we used simulations; these cause the main differences reflected in the pairs of 
values for H and K c in Tab. SI. Table SI also shows results from another study, 19 that 
employed the same kind of simulations used here and differed by obtaining X-ray data 
from oriented stacks of DMPC bilayers, from which K c was obtained directly. It also 
used the same P data, but failed to readjust the A and A values to account for the 
corrected Ka- Nevertheless, agreement is reasonable. 


Table SI: Optimal parameters found for describing the DMPC data published in Ref. 



Current 

1998a 16 

1998b 16 

2005 19 

H/ zJ 

4.11 

7.13 

4.91 

6.1 

K c /zJ 

57 ±5 

50 

80 

69 

A /Pa 

^ q 8.1 ± 0.2 

io 8 - 1 

10 81 

10 81 

A/A 

1.66 ±0.15 

1.91 

1.97 

1.91 


For completeness, the functional dependence of the individual fundamental surface 
forces for DMPC is plotted in Fig. S3. The fluctuation force becomes the dominant 
repulsive force when d\y exceeds 9 A, intermediate between the values of the Ld and 
Lo phases in Fig. 7, suggesting that the DMPC bilayer fluctuations are intermediate 
in this regard between the more fluid Ld phase and the more ordered Lo phase in the 
studied mixture. This is consistent with the Ld phase having a high concentration of 
the more disordered unsaturated lipids and the Lo phase having longer saturated chains 
with cholesterol. 


S5 SAXS analysis 

Comparisons between full q-range SAXS analyses and experimental data are shown in 
Fig. S4. Deviations between data and fits, especially for higher q ranges, are due to 
imperfect background subtraction, as explained in the section X-ray measurements in 
the main text. 
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dw( A) 

Figure S3: Partitioning of total pressure P into contributions from hydration Ph y d, van 
der Waals P v dw > and undulations P un d for DMPC. 1 ' The large open black 
circle shows the value of the separation dw at which hydration and undulation 
pressure are equal. 
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Figure S4: Calculated scattering intensities (solid lines) from full q-range analy¬ 
ses, compared to recorded SAXS data from coexisting phases (dots) of 
DOPC/DSPC/Chol (0.42:0.37:0.21), for all recorded osmotic pressures P. 
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S6 Fluctuations of the interbilayer water spacing 



Figure S5: Probability density function p of the water spacing a at different external 
pressures P, for Ld (solid) and Lo (dashed) according to Tab. 1, obtained 
from N — 32 simulations. 
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